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. Abstract Near equilibrium, the rate of relaxation to equilibrium and the transport 

(N- properties of excitations (bogolons) in a dilute Bose-Einstein condensate (BEC) are 

determined by three collision integrals, and All three collision inte- 

■ grals conserve momentum and energy during bogolon collisions, but only con- 
serves bogolon number. Previous works have considered the contribution of only two 
collision integrals, 1^^^ and J^'^. In this work, we show that the third collision integral 

. makes a significant contribution to the bogolon number relaxation rate and needs 

■ to be retained when computing relaxation properties of the BEC. We provide values 
^ , of relaxation rates in a form that can be applied to a variety of dilute Bose-Einstein 

condensates. 
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1 Introduction 

The relaxation of a dilute gas to thermal equilibrium is governed by a kinetic equation 
that contains one or more collision integrals. Each collision integral represents the ef- 
0^ ■ feet of a certain microscopic collision process on the evolution of the distribution 

function. Taken together, the collision integrals determine the rate of relaxation to 
equilibrium. In classical gases, the relaxation is governed by the Boltzmann equation 
and is determined by a single collision integral, C^^, that represents an elastic bi- 
nary collision and therefore explicitly conserves the number of particles, momentum 
and kinetic energy of the gas |[Tl|2l[3]. Similarly, for a dilute quantum gas of bosons 
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that is degenerate (but non-condensed), the relaxation is governed by the UehUng- 
Uhlenbeck equation and is determined by a modified form of C^^ that leads to a 
Bose-Einstein distribution as the equilibrium state HISl. 

Analogous kinetic equations have also been derived for dilute boson gases in 
the presence of a Bose-Einstein condensates (BEC). In ground-breaking work, Kirk- 
patrick and Dorfman |6J, derived a mean field kinetic equation (the KD equation) 
for a dilute condensed Bose gas that describes the relaxation in terms of "collisions" 
between Bogoliubov excitations. A key feature of the KD equation is that the num- 
ber of Bogliubov excitations (bogolons) need not be conserved during the collision 
process. The KD equation contains two collision integrals, which we denote as ^^^^ 
and 5^'^, that both conserve the energy and momentum of the bogolon gas. However, 
only which describes an elastic binary collision between bogolons, conserves 
bogolon number. The collision integral ^^'^ allows two (one) bogolons to enter the 
collision and one (two) to leave, thus accounting for the non-conservation of bogolon 
number during collisions. Subsequently, several authors [7 ,8 ,9| have studied trans- 
port properties of dilute BECs using mean field kinetic equations very similar to the 
KD kinetic equation with the same two collision integrals. 

We have found that there is a third collision integral, that should be included 
in the dilute BEC mean field kinetic equation, because it gives a dominate contribu- 
tion to the relaxation of the bogolon number, and a considerable contribution to the 
relaxation of other important modes in the dilute BEC. We have previously derived 
the form of the collision integral [10| and found that it can be interpreted as 
a collision where one (three) bogolons enter and three (one) leave. The purpose of 
this paper is to describe and compute the effect of the collision integral on the 
relaxation of dilute BECs to equilibrium. 

In the all of the theories described above, the relaxation to equilibrium is caused 
by the internal collision processes between particles or excitations and is described 
by collision integrals. This is distinct from the stochastic methods 111] which de- 
scribe relaxation by contact with external reservoirs, stochastic initial conditions or 
phenomenological damping. 

We begin in Section|2l with a summary of basic concepts underlying the analysis 
of the Bogolon kinetic equation. In Section [3] we write the bogolon kinetic equa- 
tion, and in Section |4] we linearize the kinetic equation in terms of deviations from 
equilibrium. In Section |5] we describe how the parameters in the collision integrals 
are obtained from a physical system. In Section |6l we obtain relaxation rates and in 
Section|2]we make some concluding remarks. 



2 Background 

We consider a spatially uniform system of bosons of mass m that are confined to 
a rectangular box of volume V with periodic boundary conditions. We assume that 
the particles interact via a contact potential y(r,, rj) = g5-'(r,- — r^), where r,- is the 
displacement of the particle and g is the strength of the interaction. 
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The Hamiltonian for this boson gas is 



1 +k2 ,k3 +k4 "^kj '^ks , 



(1) 



h k t 

where = -j;^, creates a particle with momentum /zki and at] destroys a par- 
ticle with momentum h\i\. The creation and annihilation operators satisfy the boson 
commutation relations [flkijfl^^] = 5k,,k2 where 5ki,k2 is the product of three Kjo- 
necker delta functions, one for each component of k. The summations run over all 
single particle states for both positive and negative components of k. To simplify no- 
tation in subsequent sections, we will let ^i, aj, and a^^ d\. We will 
keep this subscript convention for all quantities which are functions of the wavevector 
k/, such as Ci = and 5i,2 = 5k, .k2- 

To correctly describe a system containing a BEC, we must allow for the expecta- 
tion values (fl^flti^) and (fl_kflk) to be non-zero at equilibrium llT2l . To accomplish 
this, we introduce two mean fields, v and A, and write the unperturbed Hamiltonian 
W and interaction Hamiltonian //' as 



(£,■ - /I + v)a]ai + y (fljal,. + a^i&i) 



(2) 



and 



ijkl 



(3) 



The mean fields v and A are determined so that secular terms are eliminated from the 
kinetic equations lfT3ll . To lowest approximation, they are 



V 



(4) 



and 



(5) 



The quantity S is a shift in the energy of the system that is chosen so (//') = in 
equilibrium at zero temperature. Since S is not an operator, it will not enter into the 
collision integrals. As discussed in Ref. fW\, the chemical potential, jj., is set equal 
to V — 4 and the average particle density n = y implicity determines the values of v 
and 4. 

The mean field unperturbed Hamiltonian (|2|i can be diagonalized by a Bogohubov 
transformation. We define the operators 



= Uid^ + Vid^i 



and 



hi — Uidi + Vid; 



(6) 



with 



Ui : 



1 Ei + A 



V2 



1 Ei + A 



V2\ Ei 



(7) 
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where 

Ei = ^{Ei + Af-A^. (8) 

The operators b\ and bj have the property that [^i ,^]] — 5ij and in terms of them, the 
mean field unperturbed Hamihonian is 

H° = Y,Eib}h + WiEi-e,-A) + E. (9) 

/ ^ i 

Therefore, bJ and are interpreted as creation and annihilation operators for bosonic 
excitations, which we will refer to as bogolons. The transformation factors m, and v, 
are undefined at zero momentum, so we must give special treatment the operators 
and ao- Following Bogoliubov, we implement the special treatment of zero momen- 
tum operators by replacing aj and ao with y/No where A^o is the number of particles 
in the zero momentum state. The quantity A'o is also interpreted as the number of 
particles in the condensate, and the condensate density can be defined as no = ^• 



3 Collision Integrals for a Condensed Bose Gas 

Though we consider a dilute BEG that is spatially uniform, the momentum (energy) 
distribution of the bogolons may be out of equilibrium. The kinetic equation gov- 
erning the relaxation of the bogolon momentum distribution to equilibrium can be 
derived from the Hamiltonians presented in the preceding section. We have derived 
the mean field kinetic equation for dilute BECs using the approach of Peletminksii 
and Yatsenko (PY) II14II15II . The derivation of the collision integrals can be found in 
Ref. [ 10|. We obtain a kinetic equation for the relaxation of bogolon momentum dis- 
tribution that is nearly identical to the KD equation, but contains the a third collision 
integrals, ^l'^', in addition to the known collision integrals ^^^^ and Taking the 
thermodynamic limit by letting V ^ °° and N ^ °° while keeping « = ^ constant, 
the kinetic equation for an infinite, uniform bogolon gas can be written 

^4r^ = ^^H^} + ^^'iy^}+^l'{:^}, (10) 

where ^(k) = {P^k) is the expected number of bogolons with momentum k. The 
bogolon collision integrals are given by (we include ^^^^ and for completeness 
and comparison) 

253(ki +k2 -k3)5(£i +E2-E^){wll^)\^i^2^i-^i^2^3) 

+ 5^(ki -k2-k3)5(£i -£2-£3)(W3;|l)^(^l-^2-^-^^2^3) 

(11) 
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{Infh . 

x53(ki +k2-k3-k4)5(£i+£2-£3-£4) 



(12) 



and 



{Infh 



iik2t/k3(ik4 



^ 53 (ki - k2 - k3 - k4)5 (£1 - £2 - £3 - £4) 

X ^2,3,4)^(^1 -^2^^ - ,yK^2^3^4) 

+5^(ki+k2 + k3-k4)5(£i+£2 + £3-£4) 



(13) 



where = 1 + ._yff. 



The weighting functions are given in terms of u\ and v, by 



^l',2,3 = MIM2M3 - M1V2M3 - V1M2M3 + M1V2V3 + V1M2V3 - V1V2V3, 



(14) 



Wj^2 3 4 — M1M2M3M4 + M1V2M3V4 + U1V2V3U4 + V1M2M3V4 + V1U2V3U4 + V1V2V3V4 (15) 

and 

/31 



1,2,3,4 - 



: M1M2M3V4 + M1M2V3M4 + M1V2M3M4 + V1V2VSU4 + V1V2M3V4 + V\U2Vj,V4. (16) 



Each of these weighting functions has specific symmetry with respect to interchanges 
of its indices that is shared by its corresponding collision integral. The collision in- 
tegrals and ^^^^ are identical those given in Ref. |6| and discussed extensively 
in Ref. |i9j. It is not difficult to show that a Bose-Einstein distribution in the bogolon 
energies, 



1 



3W - 1 



(17) 



is the long-time steady-state solution of Eq. ( fTOl l. In Ref. fW\ we prove that these col- 
lision integrals conserve average particle number even through they do not conserve 
bogolon number. 



4 Linearized Collision Operators 

The relaxation rates of the bogolon gas are defined as the rates of pure exponential 
decay for infinitesimal perturbations about equilibrium. To determine the perturba- 
tions that exhibit pure exponential decay and their associated rates, we linearize the 
collision integrals about equilibrium to obtain the linearized collision operators. Re- 
placing the collision integrals in the kinetic equation with the linearized collision 
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operators results in a linear kinetic equation that can be cast as an eigenvalue equa- 
tion. The eigenvalues of this equation give the relaxation rates and the eigenfunctions 
give the characteristic relaxation modes of the perturbations. 

We linearize the collision integrals by writing the bogolon distribution as 

^it)=^° + ^''^U:{t) (18) 

where 0,(r) represents deviations from absolute equilibrium and decays to zero as 
f — ). oo. We assume that the gas is very close to equilibrium and neglect terms beyond 
first order in 0,(f). The factor of ^"^^ multiplying 0,(f) simplifies the eventual 
form of the linearized collision operators. We also introduce a dimensionless mo- 
mentum vector c = ^\ ^ k. In terms of c the integrals become dimensionless, and 
the physical quantities can be combined into an overall rate constant with units of 
inverse time. Using s = 
rate constant takes the form 



inverse time. Using g = where a is the measured s-wave scattering length, this 



^^WMT^ (19) 
nh 

To parameterize the collision integrals in terms of dimensionless quantities, we define 
b — -jAj! as the dimensionless ratio between the mean field and the temperature. Since 

'^^^ contains A^o in its coefficient, we introduce a third dimensionless quantity, 

SO that the overall coefficient of ^^^^ is ay. 

After the collision integrals have been Unearized, we can rename the integration 
variables and use the symmetries of the weighting functions ( fT4lll5lfT6b to write the 
kinetic equation as 



dt 



£/C2^°0(C2,Of2a7A(Ci,C2)-2arB(Ci,C2)-2arB(C2,Ci) 



where 



-6A(Cl,C2)-2/;A(Cl,C2)+2gB(Ci,C2) -ec(Cl,C2)-ec(C2,Cl) 



M(ci) = — ^ / t/C2^°[2arA(Ci,C2) + arB(cl,C2) + eA(Cl,C2) 



(21) 



(22) 



and 



7:i(Ci,C2)= /fl'C35^(Ci+C2-C3)5(<fi+^2-^3)W,2,3r^3, (23) 
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7i(ci,C2) - —4 / dc35\ci-C2-C3)5i<^i-^2-'^3)iWlliy^^\ (24) 



2a(ci,C2) = j dCidcAd^{c\ +C2-C3-C4) 



(25) 



Ra{Ci,C2) =-7^ / iiC3iiC45^(Ci -C2 + C3-C4) 

^2 (26) 

X - <f2 + ^3 - A)(Wf3,2,4)2^3°^4 , 

Gb(Ci,C2) = / t/C3t/C45^(Cl +C2 + C3-C4) 

J (27) 

X 5(^1 + ^2 + - A)(W^i2.1 )^-^3''^4 > 
Gc(Cl,C2)=— ^ / iiC3t/C45^(Ci-C2-C3-C4) 

J (28) 

X5(^l - ^2 - '^S - ^4)(Wi32 3 4)'=^° -^4 • 



In these expressions, (oi = ycj + 2bcf. 

The kinetic equation (l2ll i can now be cast as a general integral equation of the 



form 



Jn_±l^M{c,)Hcut)+JdC2-^ 



M(Ci)(/)(Ci,f)+ / flfC2^/:(Ci,C2)(/)(C2,f)- (29) 

1 



This equation can be simplified by expanding the angular dependence of the function 
(^(c,f) in spherical harmonics. To also make the resulting kernel symmetric under 
interchange of its arguments, we expand </>(c,f) as 



^ 00 / 



^^^^'^ ^ /^o. ^^o. ^ E WiAc,t)Yric). (30) 

Cv/^0(c)^0(c) z=Om=-/ 

Here, Y[^{c) is a spherical harmonic and \lfi^i„{c,t) is an expansion coefficient that is 
a function of c. Substitution of this expansion into Eq. ( |29T l leads to 



fl'i2iflfI22F;"'(ci)7i:(ci ,C2)F,r (C2). 

(31) 

Since K{c\ ,02) is continuous and rotationally invariant, it is a continuous function of 
ci, C2 and Ci • £2 only, and therefore can be represented by a Legendre series of the 
form 

/:(ci,C2) =£^/(ci,c2)P/(ci -02). (32) 
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Using the spherical harmonic addition theorem, we can deduce that 

1 

^ dQ,dQ2YriCl)KiCuC2)Y;r' iC2) ^27:d,j,d,„,,„, J ^/(Ci •C2)P/(Ci •C2)/:(Ci,C2) 

(33) 

and therefore Eq. dJTT l reduces to 

' "^^'di^)'^ =M{ciWl,m{cut) + j dC2K,{cuC2)WlA^2,t) (34) 

where 



K,{ci,c2) =2nc,c2J I d{ci-C2)P,{ci-C2)K{cuC2). (35) 



-1 



We have reduced the three-dimensional integral equation ( |29] l to an uncoupled set of 
one-dimensional integral equations (l34l i indexed by / and m. 

Now let us define an eigenvalue X"^^^ and an eigenfunction yf"^{c) as the pair 
which satisfies 

oo 

KmVTA'^l) ^M{C,)WIM+ j «(C1,C2)V//%(C2). (36) 



These eigenfunctions have the orthogonaUty relation 



fl'cv/LWv/," (c)-5„v. (37) 



Note that the kernel Ki{c\,C2) does not depend on m, so the eigenvalues and eigen- 
functions are also independent of m. We shall drop the m dependence of the eigenval- 
ues A",jj and eigenfunctions y/"^ (c) where it will be understood that these quantities 
are 2/ + 1 - fold degenerate. 

Using the property ( |36] | of the eigenfunctions '/"(c), Eq. (l34l i can be solved and 
the evolution of the deviation function (j){c,t) can be written in terms of the eigenval- 
ues and eigenfunctions as 

'^^'^'^^ = /x.o,\^OM Hi Al,„e-'"^y,?ic)Yr{c), (38) 

where the expansion coefficients A"^^^ depend on the initial deviation from equilibrium 
as described by (j){c,t —0). 

The eigenvalue equations for both the Boltzmann equation |2| and the Uehling- 
Uhlenbeck equation |5 1 have the same structure as Eq. ( l34l i and consequently they 
have a maximum discrete eigenvalue to which all other discrete eigenvalues converge. 
The existence of this maximum discrete eigenvalue for the Boltzmann equation has 
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been proven 1 16] and its value was shown to be equal to the minimum value of the 
function M(c). For the BEC, we denote this maximum discrete eigenvalue as Xm- 

The fact that energy is conserved implies that there is an Z = eigenfunction 
proportional to S (c) with zero eigenvalue. Similarly, the fact that momentum is con- 
served implies that there is an / = 1 eigenfunction proportional to c with zero eigen- 
value. This allows us to deduce that 

V/;L1'(c) - cS{c) y'^0(c)^0(c) (39) 

and 

Yi=i (c) - c2 ^^0(c)^0(c). (40) 

The normalization constant of these eigenfunctions must be determined numerically. 

We can use the eigenfunctions V'f (c) to write an expansion for the bogolon dis- 
tribution in the Unear regime as 

^(c) = ^O(c) + ^':me-''''w^ic)w"ic)Yr{c), (41) 

^ y'^) n=Ql=Qm=-l 

where 

<„, = / dc^ic) [^(c) - ^«(c)] (42) 
J V/q [c) 

The expansion coefficients £^qq and are exactly zero because the modes V'o('^) 
and ij/iic) are the modes associated with energy and momentum conservation and do 
not decay. Their contributions are already contained in the equilibrium distribution 
^^{c), which is a stationary state of the kinetic equation. 

5 Equilibrium Properties of the Condensed Bose Gas 

Because the collision operators are linearized about equilibrium, the eigenvalues of 
the linearized collision operators will depend on detailed properties of the equilibrium 
state. Experimentally, the equilibrium state is specified by the particle properties m 
and a, the total density n, and the temperature T. To connect our calculation to a 
physical system, we must be able to relate the parameters b and a to these physical 
properties. We do this by using equations Q and (|5]i to calculate the condensate 
density no and the mean field strength A. The dimensionless parameters b and a, 
which determine the dimensionless eigenvalues, can then be calculated from A and 
no- 

For a low density gas, the simplest consistent approximation to equations (01 and 
Q is the "Popov approximation" IqI FTTI . This approximation has been shown ifTSl 
\19\ to reproduce the experimentally measured properties of the dilute BECs very well 
for temperatures below about 60% of Tc, the critical temperature for Bose-Einstein 
condensation. For the low densities and temperatures that we consider, the mean field 
values in the Popov approximation differ from more detailed approximations ff 20ll2TI 
by only a few percent. In what follows, we use the Popov approximation to determine 
the equilibrium properties of the Bose-condensed gas. 
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Two dimensionless quantities which are readily available for this system are the 
reduced temperature T = T /Tc and the "gas parameter" Tj = na^. These quantities 
are also quite convenient for calculation, since they allow equation (HI and (|5]l to 
be written in terms of only dimensionless quantities. In the thermodynamic limit 
discussed in Sec. [3] and using the Popov approximation, Eq. (HI can be written as 



and Eq. Q can be written as 



where tjo ~ noa^, ~ (|) ~ 2.612 and 



.= |£V"^. (44) 



F{b)^^ dx \ ^. (45) 



Once Eq. (|43] | has been solved numerically for T70 with a given value of r\ and 
T, the parameters entering the collision integrals can be calculated from Eq. (l44l) and 
the expressions 

a ^ ;r3/2^!?0^-3/2 (4^) 
J? 

and 

7= 5-77^/''^ (47) 

which are obtained from Eqs. ( fT9] l. (|20] | and the definition of Tc- 

These equations show that if two systems have the same values of na^ and T /Tc, 
their condensate fractions (no/n) and their dimensionless eigenvalues will also be the 
same. However, since 7 includes m and a explicitly, the relaxation rates will differ. In 
fact, Eq. (l47T i implies an important scaling behavior: If two different systems (called 
A and B) have the same values of na^ and T /Tc, then = RBinsajj where R 

is any characteristic relaxation rate. This enables us to predict, for example, that the 
relaxation rates of ^^Rb (a = 105 gq) at T = 7.26 nK and 1.44 x lO'^ m^^ will be 
approximately 7.26% of the relaxation rates of ^^Na (a = 55 gq) at T = 0. 1 jxK and 
n= 1.0 X 10^" m^^. In this manner, we can use results at a specific value of na^ and 
T /Tc to make predictions for any system with the same values of na^ and T /Tc as 
long as its mass and scattering length are known. 

To ensure that our choice of the parameters na^ and T /Tc coincides with exper- 
imentally accessible systems, we use a dilute atomic gas of ^^Na [22] in the lower 
hyperfine state with a = 55 aq |23 | where gq is the Bohr radius. The calculations are 
carried out at a density corresponding to a gas parameter of na^ = 2.5 x 10^^. Figure 
[U shows the variation of the values of no /n, b, a and 7 with temperature as com- 
puted using the Popov approximation. For comparison, we have also plotted points 
corresponding to na^ — 1.0 x 10^^. 
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6 Calculation of Relaxation Rates 

In this section, we explain our numerical method for obtaining eigenvalues and eigen- 
vectors from Eq. ( l34l i and give the results of the calculation. To discretize Eq. ( l34l i so 
that it can be solved numerically, we use a finite ordinate method |24|. In this method, 
a quadrature scheme is chosen for evaluating integrals on the interval < c < oo and 
this scheme is applied directly to the integral in Eq. (l34l i. If we then demand that the 
eigenvalue equation is satisfied when evaluated at each quadrature point, we obtain 
the matrix equation 

Nq 

(x,) = M(x,) v/r(x,) + ]2 WjK,{x„Xi)xi/f{xj) (48) 

j=i 

where Nq is the total number of points in this quadrature and x, and w, are the quadra- 
ture points and weights, respectively. This discretization method has the advantage of 
speed and simplicity and is quite reliable, particularly if an appropriate quadrature 
scheme is chosen. For our calculations, we chose to split the integration region into 
< C2 < 2 and 2 < C2 < 12 and use a composite seven-point gaussian quadrature 
method in each region. This quadrature scheme captures the small-scale and large- 
scale features of the kernel. 

In order to gauge the accuracy of the eigenvalues we perform the calculation sev- 
eral times, each with a larger total number of quadrature points Nq. The eigenvalues 
that we report are found by performing a regression of the eigenvalues versus 1 /Nq 
and extrapolating to 1 /Nq = 0. This gives us an estimate of the truncation error in the 
eigenvalues. For all non-zero eigenvalues, this truncation error is less than 1% of the 
eigenvalue. For the zero-eigenvalues, we found that a 95% confidence interval based 
on this truncation error always includes zero. We must note that the truncation error 
only estimates the error of using a finite size matrix, and does not address the error 
involved in the computation of matrix elements. 

To calculate the matrix elements Ki{xi,Xj), we use a Gauss-Kronrod (G7-K15) 
global adaptive quadrature scheme with a relative tolerance of 10^^ Appendix 
|A] contains an in-depth explanation of how we evaluate Ki{xi,Xj) given its definition 
in Eq. ( |35] ). 

In Fig. |2] we plot several of the dimensionless eigenvalues obtained at varying 
temperature. These eigenvalues are computed for na^ = 2.5 x 10^^. Note the presence 
of a single zero eigenvalues in the / = plot, showing energy conservation, and the 
presence of another single zero eigenvalues in the / = 1 plot, showing momentum 
conservation. Analysis of the numerically calculated eigenfunctions associated with 
these zero eigenvalues supports this interpretation. The absence of any other zero 
eigenvalues confirms the fact that bogolon number is not conserved. We also find that 
the dimensionless eigenvalues appear to converge to Am from below. 

Figure |3] shows the physical relaxation rates of the some of the lower relaxation 
modes above and below Tc. The main feature is an increase in the relaxation rates 
as the temperature is lowered below Tc, followed by decrease as the temperature 
approaches zero. This increase is caused by the appearance of the condensate, and 
the start of 1 2 and 1 f^- 3 processes described by J^'^ and J^^'. The decrease as 
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the temperature approaches zero is caused by the decreasing overall rate factor 7, 
which goes as T^. For comparison, we have plotted the relaxation rates of a boson 
gas above 7c, which are obtained from a similar analysis of the Uehling-Uhlenbeck 
equation 1 5 1 . 

Since the collision integral is not included in other works, it is of interest to 
determine its effect on the eigenvalues. To do this, we have computed the eigenvalues 
without including and compared them with the eigenvalues when ^^^^ is included 
in the calculation. The results of this comparison are seen in Fig.|4l which shows the 
percent difference in the eigenvalues if the collision integral is neglected. The 
neglect of has a significant effect on the eigenvalues. As expected, neglecting 
implies neglecting 1^3 collisions and results in decreased relaxation rates for 
all modes. 

Of special note is the influence of '^^^ on Aj. To explain this, let us note that 
above the critical point, there two zero eigenvalues with / = and these are associated 
with particle number and energy conservation. Below the critical point, only one 
I — eigenvalue is zero, which must be associated with energy conservation, since 
that can be shown explicitly for the collision operator. The smallest non-zero I — 
eigenvalue tends to zero as T — Tc and is strongly affected by ^31. Since is does 
not conserve bogolon number, we interpret the smallest non-zero / = eigenvalue as 
being strongly associated with the relaxation of bogolon number 

7 Conclusions 

We have considered a new collision integral, and its influence on the relaxation 
properties of a gas of bogolons in a dilute Bose-Einstein condensate. To quantify the 
effects of this collision integral, we have computed the characteristic relaxation rates 
of the bogolon gas when it is included and when it is neglected. The new collision 
integral has a significant influence on the relaxation rates and is as important as the 
two known collision integrals in describing the relaxation of the bogolon gas. 

Linearizing the bogolon kinetic equation for small perturbations about equilib- 
rium generates a linear integral equation which can be treated as an eigenvalue equa- 
tion for the characteristic relaxation eigenmodes and relaxation rates. We have numer- 
ically computed the eigenvalues and eigenvectors of the associated matrix equation. 
These eigenvalues and eigenvectors give us the characteristic relaxation rates of the 
bogolon gas and a mode expansion for the relaxation of the momentum distribution to 
equilibrium. We find that the presence of four zero eigenvalues supports the fact that 
energy and momentum are conserved during during every bogolon collision, while 
the absence of a fifth zero eigenvalue supports the fact that bogolon number is not 
conserved. 

The characteristic relaxation rates can be related to the transport coefficients. In 
general, the expressions for the transport coefficients involve a sum over the inverses 
of the relaxation rates, so no single relaxation rate determines a transport coefficient. 
From the angular dependence of the eigenfunctions we can make some conclusions 
based on analogies to the classical Boltzmann Equation |2|. Since shear viscosity is 
related to momentum current correlations, which have a dyadic character, we expect 
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that the I = 2 eigenvalues will determine the shear viscosity. Similarly, we expect that 
the / = 1 eigenvalues will determine the thermal conductivity, since it is related to en- 
ergy current correlations, which have a vector character We expect that the remaining 
viscosities of the BEC will depend on the I = eigenvalues. 
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A The Form of the Kernels used in Computation 

To calculate the values of the kernel Ki (ci , C2 ) in Eq. i35\ . we split the calculation into six parts, one for 
each of the individual kernels in Eqs. )25I - I28t . In this appendix, we show how to obtain an expression that 
is well-suited to numerical quadrature for the kernel 2^(ci ,C2). Similar procedures can be used with the 
other five kernels as well. We begin with the definition 



J (49) 

53(c, +C2 -C3 -C4)5(f?, +<f2 -'^i -<«l)(»'i^i3,4)^^3'^4- 
First let us perform the integration over C4, 



ei(c,,C2) =271 1 d{ci ■C2)P,(C, -£2) I ^/C35(^, +^?2 - <^3 - <g»)K2,3,4)'^3'^4 ■ 



(50) 



The resulting integrand only depends on the magnitude of C4 , where C4 = | C| + C2 — C3 1 . We now perform 
the integration over C3 in spherical coordinates, with the z-axis oriented parallel to Ci +C2, 



(51) 



^0 

1 ■^A ■ 



Q'a{<:uC2)=27C J d{c, ■C2)P/(ci -£2) J 4dc3 
-1 
1 2iz 

xj dzij d035(<?l+^2-4-'«t)(Wf2,3,4)-^3' 
-1 

Now notice that our choice of spherical coordinates for the C3 integration allows us to write C4 as C4 = 
y^|ci + C2p + Cj — 2c3 |ci +C2IZ3. We can change variables to write the 23 integration as an integration 
over C4 with dc^ = — ^^^J|jj^(iz3, 



QA{ci,c2)='i7c^ I d{crC2)Pi{crC2) I clda 
-1 b 
IC1+C2I+C3 

X / dC4 , g(<gl+^2-^3-A)(H^p2 3 4)'-y3°-^4° 

./ C3|Ci +C2I 
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Since the integrand now only depends on ci, and Ci ■ £2, let us define ca = \J c\ + 0^ + 2ciC2(ci £2) 
and use ca as a change of variables for the £1 ■ £2 integration. This results in 

QfA{cuC2)=An^ TdcA^P,('^-^) 
J ciC2 \ 2c|C2 / 

(53) 

X / cjdc^ f dc^ — S(Si+S2-Si-<S^) {wH 3 4)^^° ^4 ■ 

J J CiCA ' ' ' 

H-Cil 

To handle the integration over ca, we write the integration limits in terms of Heaviside theta functions, 



Q'a(<'-1,C2)=471^-^I CidCij C^dC4 I dCAPl {—^^1^^ S{^l+^2-<^3-£4) 





x(W'L2.3.4)^^3°-^4e(kA -C3I < C4 < Ca + C3)6{\ci - nl < Q < Cl +q) 
and notice that 

e(|cA -C3I < C4 < CA + C3) = 9(\C3 - C4\ < CA < C3 + C4). (55) 
This allows us to move the ca integration through all of the others and write 

Q'a(ci,C2)=47Z^— I C3dC3 I C4dC4S(Sl+S2-<Si-<g4)(Wll34f.^l^Wl2iA, (56) 
C\C2 J .1 ' ' ' 



(54) 



where 



""1,2.3,4 - 



minkl+Q.O+Ql , j 2\ 



raax||r|-C2|,k3-C4|| 



The final delta function of energy can now be handled in several ways, but each of them will lead to a 
well-behaved integrand. Still, we can make a few observations that will help the quadrature go faster 

First, notice that the range of integration on C3 must satisfy S-i < S\+S2. This corresponds to C3 < 



' ^ {£'1 + ^2)^ + — b. In fact, the symmetry between C3 and C4 shows that the whole integral is equal 



to twice the integral from < C3 < c/, where '■I' ~ ]J \J ( ^^t^) +lP- — b. Also, since it is symmetric in 

ci and Q, we can assume that c\ > q in Eq. (57) and swap ci and C2 if £'2 > ■ 

Analysis of the function n'[ 234 shows that under the constraint that ^1 + (02 = (?3 + A, 2 3 4 = 
2min[ci,C2,C3,C4]. This is only generally true for and not the other kernels. Furthermore, W[234 
for / > 1 can always be written as the product of Wj j 3 4 and another finite function. The integrand will 
therefore always have a discontinuity of its derivative when C3 = C2, and the integration region should be 
spHt at this point. 

With all of these considerations, we finally write the best form as 



Q'AicuC2) = — 



j dci + j dc3 



^(Wa3,4)'^3''^4V,2,3,4 (58) 
cj + b 



where C4 takes the value that makes tSj = tf 1 + (^2 ^ ■ Three of the other kernels (Qg, Qc, Ra) can be 
similarly reduced to a single integration with a well-behaved integrand, while the other two (Ta, Tg) can 
be reduced to explicit functions of c\ and C2. 

A final point concerning the kernels involves the function Rjj(ci,C2). Though this function is un- 
defined when c\ = C2, integrals over the entire kernel Ki{c[,C2) still converge. In the numerical method 
of Sec. (6), we use the fact that ^^'^^ acting alone conserves bogolon number to determine the values of 
R°a(^>.ci). 
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Fig. 1 Variation of the condensate fraction, dimensionless parameters h and a and the overall rate coeffi- 
cient 7 for -^Na as the temperature is increased. All quantities are computed in the Popov approximation. 
Circles represent a gas parameter of na^ = 2.5 x 10^* and squares represent no' = 1.0 x 10^^. The solid 
line is the condensate fraction of an ideal Bose gas. 
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Fig. 2 Dimensionless eigenvalues X" as a function of reduced temperature T jTc for / = (a), / = 1 (b) 
and 1 = 2 (c). Solid circles indicate the maximum discrete eigenvalue Xm- AU quantities are computed in 
the Popov approximation at nc? = 2.5 x 10~*. 
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Fig. 3 Relaxation rates ■yX" versus temperature for several of the lowest eigenvalues. All quantities are 
computed in the Popov approximation for ^'Na at no' = 2.5 x 10~*. For comparison, relaxation rates of 
the Uehling-Uhlenbeck equation are plotted for T > Tc. 
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Fig. 4 Percentage change in several of the lowest non-zero relaxation rates if the effects of are ne- 
glected. Relaxation rates is significantly slower when is neglected, especially for Xq , which is the 
most significant eigenvalue in the relaxation of bogolon number. 



